APS/123-QED 

Directionality is an inherent property of biochemical networks 

Feng YangJ*| Feng Qi, and Daniel A. Bearcl 
Biotechnology and Bioengineering Center, Department of Physiology, 
Medical College of Wisconsin, Milwaukee, Wisconsin 53226, USA 
! (Dated: February 9, 2008) 

^ ■ Abstract 

Thermodynamic constraints on reactions directions are inherent in the structure of a given 

biochemical network. However, concrete procedures for determining feasible reaction directions for 

large-scale metabolic networks are not well established. This work introduces a systematic approach 

to compute reaction directions, which are constrained by mass balance and thermodynamics, for 

genome-scale networks. In addition, it is shown that the nonconvex solution space constrained by 

physico chemical constraints can be approximated by a set of linearized subspaces in which mass 

and thermodynamic balance are guaranteed. The developed methodology can be used to ab initio 

£\\ , predict reaction directions of genome-scale networks based solely on the network stoichoimetry. 
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Constraint-based approaches to analyzing biochemical networks, such as flux balance 



analysis (FBA) and energy balance ana 
system analysis of metabolic networks 



ysis (EBA), have found widespread applications in 
l|. A successful set of procedures has been estab- 



lished for determining metabolic reactions which are present in a genome-scale system 
However, the procedure for determining feasible reaction directions is less concrete and ar- 
bitrary to some degree . Yet the information on reaction directions is crucial to ensure 
that predicted reaction fluxes are physically feasible and obtained results are physiologically 
reasonable 

Allowable reaction directions in a biochemical network are constrained by mass- 
conservation constraints and thermodynamic constraints, both arising from the stoichiome- 

nnnu 

try of a given network [5|, 10, [8|, 19[ . While the linear mass-conservation constraint is rigorously 
applied in many applications, the nonlinear thermodynamic constraint is not typically con- 
sidered due to its implicit nonlinearity and NP-completeness of the associated problems. 
Consequently, the flux distribution determined solely by the principle of mass-conservation, 
may not be physically feasible. 

The nonlinear thermodynamic constraint is based on the second law of thermodynamics 
which requires that each internal reaction with non-zero flux must dissipate energy [71]. For 
a set of reaction directions to be thermodynamically feasible (T-feasible), there must exist 
a thermodynamic driving force for the reaction directions on a given network. A robust 
algorithm to implement this constraint, which is based on an NP-complete computation of 
entire elementary modes for a network, has been previously developed [10]. However, an 
appropriate method for genome-scale networks has not been developed until now. Here, 
we introduce a new methodology for systematically determining feasible reaction directions 
for large-scale networks and apply it to a biochemical network model of the mammalian 
cardiomyocyte [ill ]. 

To illustrate mass-balanced and thermodynamically-defined determination of feasible re- 
action directions, an example network illustrated in Fig. 1 is used. In this network, transport 
fluxes Ji, J12, and J13 are assumed to have net flux in the direction indicated by the arrows 
in the figure (defined as positive). In general, we assume that some subset of reaction fluxes 
(typically transport fluxes) have associated with them prior assigned directions. A set of 
these predefined constraints is denoted as For example of Fig. 1, the constraints on 
boundary fluxes can be expressed sign(Ji, Ji 2 , J13) = {+, +, +}. 
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FIG. 1: A example network composed of 13 reactions (J1-13) and 9 species (A through I). J\, 
J12 and J13 are transport (boundary) fluxes that transport mass into and out of the system; J2 
through Jn are internal fluxes, as shown in the dashed box. 

Flux directions constrained by mass-balance only For a flux vector J to be feasible, 
steady- state mass balance requires 

Sij ■ Jj = 0, and J obeys ^ (1) 
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where is the stoichiometric coefficient of metabolite i and reaction j. Mass-irreversible 
(MI) fluxes are those fluxes, the signs of which are specified by considering the mass-balance 
and boundary constraints of Eq. 1 alone. For example of Fig. 1, J 5 is mass- irreversible, 
because there exists no feasible flux vector J that satisfies Eq. 1 for which J 5 < 0. In 
fact, here we say flux J 5 is strictly mass-irreversible (sMI) because J5 > and the equality 
J5 = is not feasible for finite boundary flux values. Such fluxes represent the backbone 
of the network; without them, the network would not function with any nonzero fluxes. 
In addition, J7 is characterized as mass-infeasible (MN) because no flux is possible for 
this network given the boundary conditions and the network stoichiometry. The remaining 
internal reactions in Fig. 1 (J2-4, Jq, and Jg-n) are denoted as mass-reversible (MR). 
Flux directions constrained by thermodynamic constraints The nonlinear ther- 

— * 

modynamic constraint requires that there exists a driving force jl for a given flux vector J 
such that 

Ayt/j — fa ■ S^, Jj ■ Afij < 0, 

i 

and Jj = if and only if A/ij = 0, (2) 

where i and j are the indexes for a particular metabolite and internal reaction, respectively. 
Thermodynamically irreversible (TI) fluxes are those identified as mass-reversible by Eq. 1 
but found to be irreversible when T-constraints (Eq. 2) is imposed. In Fig. 1, J 6 , J 8 and 
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Jio are TI fluxes, with their feasible net fluxes in the positive directions illustrated in the 
figure. There exists no T-feasible flux patterns for which these reactions operate with net 
fluxes in the negative direction. In addition, reactions J2, J3 and J4 are characterized as 
T-infeasible (TN) because no net flux is possible in any of flux directions; reactions Jg and 
Jn can proceed in both directions while still satisfying mass-balance and T-constraints, and 
thus they are denoted as T-reversible (TR) reactions. An algorithm for determining TI, 
TN, and TR fluxes in a given network is presented in the Supplemental Material. 

In characterizing the impacts of applying T-constraints, we observe that a mass-reversible 
(MR) reaction is potentially reducible by T-constraints, either reduced as T-irreversible (TI) 
or reduced as T-infeasible (TN) . In addition, a mass-irreversible (MI) reaction is also 
reducible if it is later characterized as T-infeasible, though it may not be necessarily true. 
Furthermore, the constraints of Eq. 2 can be equivalently represented by the orthogonality 
of a flux vector J against the sign patterns of the exhaustive set of infeasible reaction 

Q 

cycles [7jj. Here, a reaction cycle is a set of reactions, the summation of which results in 
cancellation of all participating metabolites. There are 4 reaction cycles imbedded in the 
example network (Fig. 1), for example, three fluxes of J2, J3, and J4 make up a cycle: 
J 2 : E ^ F; J 3 : F ^ G; J 4 : G ^ E, and thus J 2 + J 3 + J 4 : = 0. 

Flux directions that are sufficient for T-feasibility The overall goal of ab initio 
prediction is to generate a set of irreversibility constraints that guarantees the T-feasibility. 
However, the constraints defined so far (M I and TI) are merely necessary, but not sufficient 
to ensure the T-feasibility in general [10]. Consequently, the solution space defined by MI 
and TI is not linearly constrained and still a nonconvex space. For instance, to generate a 
linearly constrained, T-feasible solution space for the example network, one of the following 
additional linear constraints has to be included: I) If Jg > 0, then Jn is feasible in any 
directions; II) If Jn < 0, then J 9 is feasible in any directions; III) If J 9 < and Jn < 0; IV) 
If Jg > and Jn > 0, etc. The resulting convex spaces (iEBA), as well as traditional flux 
balance space (FBA) and nonconvex T-feasible space(fEBA), are illustrated in Fig. 2. 

In general, in addition to those necessary conditions we need to determine a set of extra 
irreversibility constraints (EI) to guarantee the T-feasibility. The resulting linearized so- 
lution spaces (iEBA) defined by MI, TI, and EI are subsets of the nonconvex T-feasible 
space (fEBA). However, a complete enumeration of these subspaces is not practicable 
(maximally 21 n , where n is the number of T-reversible fluxes in a given network), especially 
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FIG. 2: Illustration of the solution space constrained by mass-balance only (FBA), the space 
constrained by additional nonlinear T-constraints (fEBA), and the linearlized T- feasible spaces 
(lEBA, such as the solution spaces of I, II, III, IV, etc). The space defined by traditional FBA 
may not be physically feasible; the space of fEBA, in which both mass-balance and energy-balance 
are guaranteed, represents a better approximation of network behaviors than using FBA only. 

for large-scale networks. To approximate the nonconvex fEBA space, however, it is possible 
to enumerate a limited set of subspaces (iEBA) in which the number of reversible reactions 
is maximized. The subspaces obtained represent the largest subsets of fEBA, such as the 
spaces of I and II in Fig. 2. Therefore, the network behaviors could be maximally captured 
by studying these linearized subspaces. For those subspaces that are heavily constrained, 
for instance, the spaces of III, IV and many others, it is assumed that they represent the 
non-typical network behaviors, or the time staying in these states is highly limited. In prac- 
tice, determining which set of \EBA subspaces for particular studies may require further 
constraints or evidences from experimental observations. 

In summary, for a given network with associated prior reaction directions, to generate a 
convex, T-feasible solution space, reaction directions can be categorized as {MN, sMI, MI, 
TN, TI, TR, and EI}. A diagram which shows their relationship are illustrated in Fig. 3. 
In particular, when all reaction directions other than those of predefined boundary fluxes 
are assumed reversible, a few theorems have been observed in the above analysis. 

Theorem 1 Mass-irreversible reactions are not reducible by thermodynamic 
constraints. Suppose a MI reaction that connects species A to B, A — > B, is reduced 
as T-infeasible by a reaction loop which connects species B to A along L, B A. 
As a result, there must be a pathway that connects A to B along a reversed loop of L, 
A ^~ L ^ B, by the definition of mass balance. The existence of such a reversed loop thus 
disables the reducibility of the loop L, and there is no cycle existing between species A and 
B. This theorem implies a mass-irreversible reaction does not participate in any reaction 
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FIG. 3: Illustration of the categories of reaction directions (shaded in gray) when both mass- 
conservation and T-constraints are imposed. 

cycle; if involved, then it must be characterized as mass-reversible by Eq. 1. In other words, 
MI n TI = 0, where MI and TI represent the sets of mass- irreversible and T-irreversible 
reactions, respectively. 

Theorem 2 Mass-reversible reactions can be reduced as T-irreversible or T- 
infeasible. Supposing a reversible flux Jj, A t> B, involved in a loop L, then all reactions 
in this loop L must be reversible because of the connectivity between A and B along L. To 
block this reaction cycle, any reversible reactions can be chosen and reduced as irreversible. 
Therefore, the loop constraint is active with respect to the reversible flux Jj, which can be 
reduced as T-infeasible (such as J2-4), or reduced as T-irreversible (such as J6,8,io in the 
example network). 

Theorem 3 A reaction cycle in which all reactions are irreversible is impossible. 

Suppose there is a reaction cycle which contains an irreversible reaction A — > B, and a set 
of irreversible reactions along loop L, then the reactions along loop L must be reversible by 
mass-balance. Therefore, a fully directed cycle is impossible and any cycles can always be 
eliminated by assigning directions to some T-reversible reactions in the network. 

The developed methodology has been applied for the metabolic network of mammalian 



cardiomyocyte [Uj which accounts for 240 metabolites and 257 reactions. While the stoi- 
chiometric structure of this model has been adopted, the flux directions previously assigned 
are typically not considered. Instead, we define a network in which the directions of all 
boundary fluxes and two internal fluxes (ATP hydrolysis and a transport flux of 12 DGRt2), 
are assigned based on the published model [111 ] . but all other reactions are assumed to be re- 
versible in default. In this network model, three transport fluxes are constrained to operate 
in the direction of net uptake of substrates (glucose, oleate and oxygen) and 28 transport 
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FIG. 4: Illustration of a thermodynamically feasible network of the cardiomyocyte network. The 
nodes and edges represent metabolites and reactions, respectively. In particular, the red, blue and 
yellow nodes represent the sources, the sinks and the intermediate metabolites of the network, 
respectively. Red edges and blue edges represent the strictly mass-irreversible reactions, and violet 
red edges and navy blue edges represent mass-irreversible reactions. In particular, reactions in red 
or violet red proceed in the same direction as the original definition of the network U| ; reactions 



in blue and navy blue proceed in a reversed direction. An infeasible reaction cycle is shown at the 
right top corner of Fig. 4. For simplicity, the currency of the biochemical networks, such as ATP, 
ADP, Pi, NADH, NAD, H+, H 2 and C0 2 , have not been shown. 

fluxes (including citrate, acetate and 26 additional output fluxes) are constrained to operate 
in the direction of net production of reactants. A total of 33 irreversibility constraints are 
specified and the directions of the remaining 217 internal reactions are computed. 

Given finite uptakes of glucose, oleate and oxygen, mass-balance analysis identifies 18 
sMI internal fluxes (see red and blue fluxes in Fig. 4), which represent a set of reactions 
essential for this network. It also shows that only 14 products (of 32 desired metabolites) 
could be possibly produced from the sources (i.e. glucose, oleate, and 2 ), with the aids 
of three reversible transport fluxes of C0 2 , H + , and H 2 0. In addition, 70 mass-infeasible 
fluxes, 73 mass-irreversible fluxes, and 74 mass-reversible fluxes are found. Furthermore, 
the results of intensively sampled T-feasible networks (559,518,687 samples) further indicate 
that among 74 mass-reversible fluxes, three fluxes (CITRtm, ORNt3m, ORNtAm) are 
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T-infeasible, because these three reactions 

CITRtm : citrL[m] ^ citrL[c], 

ORNt3m : h[c] + orn[m] ^ h[m] + orn[c], 

ORNtAm : citrL[c] + h[c] + orn[m] ^ cztrL[m] + /i[m] 

+orn[c] (3) 



develop a reaction cycle as C IT 'Rtm + O RNtAm — RNt3m = 0. It is also found that h[c] is 
the only metabolite that connects this reaction cycle to the remaining part of the network, 
which reminds of the reaction cycle J% — > J3 — > J4 — ► J2 in the example network. Thus, 
neither direction is feasible for all these three reactions. 

The solution space defined by MI U TI needs to be further refined since 1099 cycles 
are still embedded in this network. However, by specifying only 14 reversible fluxes (out of 
71 T-reversible fluxes) be irreversible, it is possible to eliminate all those cycles. A final T- 
feasible network for mammalian cardiomyocyte has been illustrated in Fig. 4, which consists 
of 112 infeasible reactions, 101 irreversible reactions and 44 reversible reactions. However, 
such a linearly constrained subspace is not unique, as indicated in the example network. 
To further explore the nonlinear T-feasible space, a large number of subspaces (see Fig. 
5) have been generated by assigning directions to those T-reversible fluxes. Importantly, 
each irreversibility assignment is generated in such a way that the maximum number of 
reaction cycles can be eliminated. Therefore, the generated set of linear constraints defines 
a feasible solution space which requires a minimal number of irreversibility constraints, and 
thus represents the largest linear subspace in Fig. 2. 

Surprisingly, the computed T-feasible irreversibility constraints show significant differ- 
ence when compared to the originally assigned directions, by which the space defined is not 
thermodynamically feasible (7 reaction cycles are included). To make it physically feasible, 
three mass-reversible fluxes (CITRtm, ORNt3m, ORNtAm) have to be eliminated, and at 
least one additional direction assignment, either specifying AKGMALtm (previously as- 
signed as reversible) be positive or specifying ASPTAm be negative, has to be incorporated 



into the original definition of the network. Althoug 
voke significant impacts on the overall results from 



1 these additional constraints do not in- 



Ill], they do change the flux distribution 



with the given boundary fluxes. While the revised network represents a special subset of 
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FIG. 5: Demonstration of reductions in the number of reaction cycles with extra irreversibility 
constraints (EI). It represents the results of 10 5 independent simulations for generating T- feasible 
solution spaces. Bars represent the standard deviation of the number of cycles that are eliminated 
at each particular assignment. 

TABLE I: Comparison of a T-feasible network with high flexibility to the original network as well 
as its corresponding mass-balanced and thermodynamically balanced networks. 



PI PI+FBA PI+FBA+EBA EBA 

Infeasible fluxes 149 149 112 

Positive 6 fluxes 156 60 61 42 

Negative 6 fluxes 3 21 21 59 

Reversible fluxes 98 27 26 44 



Q A T-feasible network corresponds to the network shown in Fig. 4. 

fc The positive reaction direction is defined based on the left-to-right reaction direction in the original model 

HI 

definitions [111 ]. 

the T-feasible solution space, it is significantly different from the network which shows the 
maximum level of flexibility (i.e. containing the maximum number of reversible reactions). 
Table 1 shows the difference between an example of the largest T-feasible subspaces, and 
the original definition of the network (PI) and its mass-balanced network (PI+FBA) and 
thermodynamically balanced network (PI+FBA+EBA). The significant difference between 
the ab initio prediction and the previously assigned directions, can be explained by the fact 
that mass- conservation and thermodynamic constraints are by far not the only constraints 
that a given network must satisfy. Those nonidentified may come from gene regulatory 
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constraints, topological constraints, environmental constraints [Tj] , and constraints that are 
accumulated from evolution and natural selection. By systematically studying the differ- 
ence between the results from ab initio predictions and experimental observations, we may 
find some interesting knowledge "gap", leading to more accurate predictions on the flux 
distributions. 

In summary, a biophysical constraint is applied to genome-scale biochemical net- 
works for computing flux directions and determining a linearly constrained, mass- and 
thermodynamically-balanced solution space. Given a reaction network and a set of prior 
constraints on certain flux directions (such as known transport flux directions) the algo- 
rithm determines a minimal set of irreversibility constraints (reactions directions) that is 
sufficient to guarantee T-feasibility. The developed methods could be used as a tool to sys- 
tematically design a biochemical network, greatly facilitate the further flux balance analysis 
and energy balance analysis, improve their predictions of the flux distribution, and help 
capturing the physically meaningful behavior of the network studied. Numerically, a linear 
solution space will make related optimizations more reliable and efficient. 
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Constraint-based approaches to analyzing biochemical 
networks, such as flux balance analysis (FBA) and en- 
ergy balance analysis (EBA) , have found widespread ap- 
plications in system analysis of metabolic networks [? ] . 
A successful set of procedures has been established for 
determining metabolic reactions which are present in a 
genome-scale system [? ? ]. However, the procedure for 
determining feasible reaction directions is less concrete 
and arbitrary to some degree [? ? ] . Yet the information 
on reaction directions is crucial to ensure that predicted 
reaction fluxes are physically feasible and obtained re- 
sults are physiologically reasonable [? ? ] . 

Allowable reaction directions in a biochemical net- 
work are constrained by mass-conservation constraints 
and thermodynamic constraints, both arising from the 
stoichiometry of a given network [????]. While 
the linear mass-conservation constraint is rigorously ap- 
plied in many applications, the nonlinear thermodynamic 
constraint is not typically considered due to its implicit 
nonlinearity and NP-completeness of the associated prob- 
lems. Consequently, the flux distribution determined 
solely by the principle of mass-conservation, may not be 
physically feasible. 

The nonlinear thermodynamic constraint is based on 
the second law of thermodynamics which requires that 
each internal reaction with non-zero flux must dissipate 
energy [? ]. For a set of reaction directions to be ther- 
modynamically feasible (T- feasible), there must exist a 
thermodynamic driving force for the reaction directions 
on a given network. A robust algorithm to implement 
this constraint, which is based on an NP-complete com- 
putation of entire elementary modes for a network, has 
been previously developed [? ]. However, an appropriate 
method for genome-scale networks has not been devel- 
oped until now. Here, we introduce a new methodol- 
ogy for systematically determining feasible reaction di- 
rections for large-scale networks and apply it to a bio- 




FIG. 1: A example network composed of 13 reactions (Ji_i3) 
and 9 species (A through I). Ji, J12 and J13 are transport 
(boundary) fluxes that transport mass into and out of the 
system; J2 through Jn are internal fluxes, as shown in the 
dashed box. 



chemical network model of the mammalian cardiomy- 
ocyte [? ]. 

To illustrate mass-balanced and thermodynamically- 
defined determination of feasible reaction directions, an 
example network illustrated in Fig. 1 is used. In this 
network, transport fluxes J\, J12, and J13 are assumed 
to have net flux in the direction indicated by the arrows 
in the figure (defined as positive) . In general, we assume 
that some subset of reaction fluxes (typically transport 
fluxes) have associated with them prior assigned direc- 
tions. A set of these predefined constraints is denoted as 
<5 . For example of Fig. 1, the constraints on boundary 
fluxes can be expressed ^: sign(Ji, J i2 , J13) = {+, +, 
+}■ 

Flux directions constrained by mass-balance 
only For a flux vector J to be feasible, steady-state mass 
balance requires 



Sij ■ Jj = 0, and J obeys \1/ 



(1) 



where Sy is the stoichiometric coefficient of metabolite i 
and reaction j. Mass-irreversible (MI) fluxes are those 
fluxes, the signs of which are specified by considering the 
mass-balance and boundary constraints of Eq. 1 alone. 

For example of Fig. 1, J5 is mass- irreversible, because 

there exists no feasible flux vector J that satisfies Eq. 1 

'Electronic address: dbeard@mcw.edu. However, Daniel A. Beard doeslMtW^ifihwi'iEji <h^ovkmta^ailimtimCa^^i4hmiet^ iflarflsHri^. 
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mass-irreversible (sMI) because J5 > and the equality 
J5 = is not feasible for finite boundary flux values. Such 
fluxes represent the backbone of the network; without 
them, the network would not function with any nonzero 
fluxes. In addition, J7 is characterized as mass-infeasible 
(MN) because no flux is possible for this network given 
the boundary conditions and the network stoichiomctry. 
The remaining internal reactions in Fig. 1 (J2-4, J§, and 
J8-11) are denoted as mass-reversible (MR). 

Flux directions constrained by thermodynamic 
constraints The nonlinear thermodynamic constraint 
requires that there exists a driving force jl for a given 
flux vector J such that 

A/ij = J^Mj ' Sij, Jj ■ Afij < 0, 

i 

and Jj = if and only if A^ij = 0, (2) 

where i and j are the indexes for a particular metabo- 
lite and internal reaction, respectively. Thermodynami- 
cally irreversible (TI) fluxes are those identified as mass- 
reversible by Eq. 1 but found to be irreversible when 
T-constraints (Eq. 2) is imposed. In Fig. 1, J 6 , Js and 
J10 are TI fluxes, with their feasible net fluxes in the pos- 
itive directions illustrated in the figure. There exists no 
T-feasible flux patterns for which these reactions operate 
with net fluxes in the negative direction. In addition, 
reactions J2, J3 and J 4 are characterized as T-infeasible 
(TN) because no net flux is possible in any of flux direc- 
tions; reactions Jg and Jh can proceed in both directions 
while still satisfying mass-balance and T-constraints, and 
thus they are denoted as T-reversible (TR) reactions. An 
algorithm for determining TI, TN, and TR fluxes in a 
given network is presented in the Supplemental Material. 

In characterizing the impacts of applying T- 
constraints, we observe that a mass-reversible (MR) re- 
action is potentially reducible by T-constraints, either re- 
duced as T-irreversible (TI) or reduced as T-infeasible 
(TN) . In addition, a mass-irreversible (MI) reaction is 
also reducible if it is later characterized as T-infeasible, 
though it may not be necessarily true. Furthermore, 
the constraints of Eq. 2 can be equivalently represented 
by the orthogonality of a flux vector J against the sign 
patterns of the exhaustive set of infeasible reaction cy- 
cles [? ]. Here, a reaction cycle is a set of reactions, 
the summation of which results in cancellation of all 
participating metabolites. There are 4 reaction cycles 
imbedded in the example network (Fig. 1), for exam- 
ple, three fluxes of J2, J3, and J 4 make up a cycle: 
J 2 : E ^ F; J 3 : F ^ G; J 4 : G ^ E, and thus 

J2 + J3 + J4 ■ = 0. 

Flux directions that are sufficient for T- 
feasibility The overall goal of ab initio prediction is to 
generate a set of irreversibility constraints that guaran- 
tees the T-fcasibility. However, the constraints defined so 
far (MI and TI) are merely necessary, but not sufficient 
to ensure the T- feasibility in general [? ]. Consequently, 
the solution space defined by MI and TI is not linearly 
constrained and still a nonconvex space. For instance, to 
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FIG. 2: Illustration of the solution space constrained by mass- 
balance only ( FBA) , the space constrained by additional non- 
linear T-constraints (fEBA), and the linearlized T-feasible 
spaces (iEBA, such as the solution spaces of I, II, III, IV, etc). 
The space defined by traditional FBA may not be physically 
feasible; the space of fEBA, in which both mass-balance and 
energy-balance are guaranteed, represents a better approxi- 
mation of network behaviors than using FBA only. 



generate a linearly constrained, T-feasible solution space 
for the example network, one of the following additional 
linear constraints has to be included: I) If Jg > 0, then 
Jn is feasible in any directions; II) If Jn < 0, then Jg is 
feasible in any directions; III) If Jg < and Jn < 0; IV) 
If Jg > and Jn > 0, etc. The resulting convex spaces 
(iEBA), as well as traditional flux balance space (FBA) 
and nonconvex T-feasible space(fEBA), are illustrated 
in Fig. 2. 

In general, in addition to those necessary conditions 
we need to determine a set of extra irreversibility con- 
straints (EI) to guarantee the T-feasibility. The result- 
ing linearized solution spaces (iEBA) defined by MI, 
TI, and EI are subsets of the nonconvex T-feasible space 
(fEBA). However, a complete enumeration of these sub- 
spaces is not practicable (maximally 21™, where n is the 
number of T-reversible fluxes in a given network), espe- 
cially for large-scale networks. To approximate the non- 
convex fEBA space, however, it is possible to enumerate 
a limited set of subspaces (iEBA) in which the number 
of reversible reactions is maximized. The subspaces ob- 
tained represent the largest subsets of fEBA, such as the 
spaces of I and II in Fig. 2. Therefore, the network be- 
haviors could be maximally captured by studying these 
linearized subspaces. For those subspaces that are heav- 
ily constrained, for instance, the spaces of III, IV and 
many others, it is assumed that they represent the non- 
typical network behaviors, or the time staying in these 
states is highly limited. In practice, determining which 
set of iEBA subspaces for particular studies may require 
further constraints or evidences from experimental obser- 
vations. 

In summary, for a given network with associated prior 
reaction directions, to generate a convex, T-feasible so- 
lution space, reaction directions can be categorized as 
{MN, sMI, MI, TN, TI, TR, and EI}. A diagram which 
shows their relationship are illustrated in Fig. 3. In par- 
ticular, when all reaction directions other than those of 
predefined boundary fluxes are assumed reversible, a few 
theorems have been observed in the above analysis. 
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FIG. 3: Illustration of the categories of reaction direc- 
tions (shaded in gray) when both mass-conservation and T- 
constraints are imposed. 



Theorem 1 Mass-irreversible reactions are not 
reducible by thermodynamic constraints. Suppose 
a MI reaction that connects species A to B, A — > B, is 
reduced as T-infeasible by a reaction loop which connects 
species B to A along L, B -^ L ^ A. As a result, there 
must be a pathway that connects A to B along a reversed 
loop of L, A -^>~ L ^r B, by the definition of mass balance. 
The existence of such a reversed loop thus disables the 
rcducibility of the loop L, and there is no cycle existing 
between species A and B. This theorem implies a mass- 
irreversible reaction does not participate in any reaction 
cycle; if involved, then it must be characterized as mass- 
reversible by Eq. 1. In other words, MI n TI = 0, 
where MI and TI represent the sets of mass-irreversible 
and T- irreversible reactions, respectively. 

Theorem 2 Mass-reversible reactions can be re- 
duced as T-irreversible or T-infeasible. Supposing 
a reversible flux J t , A ^5 B, involved in a loop L, then 
all reactions in this loop L must be reversible because of 
the connectivity between A and B along L. To block this 
reaction cycle, any reversible reactions can be chosen and 
reduced as irreversible. Therefore, the loop constraint is 
active with respect to the reversible flux Jj, which can 
be reduced as T-infeasible (such as J2-4), or reduced as 
T-irreversible (such as J6,8,io in the example network). 

Theorem 3 A reaction cycle in which all reac- 
tions are irreversible is impossible. Suppose there 
is a reaction cycle which contains an irreversible reaction 
A — > B, and a set of irreversible reactions along loop 
L, then the reactions along loop L must be reversible 
by mass-balance. Therefore, a fully directed cycle is im- 
possible and any cycles can always be eliminated by as- 
signing directions to some T-reversible reactions in the 
network. 

The developed methodology has been applied for the 
metabolic network of mammalian cardiomyocyte [? | 
which accounts for 240 metabolites and 257 reactions. 
While the stoichiometric structure of this model has been 
adopted, the flux directions previously assigned are typ- 
ically not considered. Instead, we define a network in 
which the directions of all boundary fluxes and two in- 
ternal fluxes (ATP hydrolysis and a transport flux of 
l2DGRt2), are assigned based on the published model 
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FIG. 4: Illustration of a thermodynamicalfy feasible network 
of the cardiomyocyte network. The nodes and edges represent 
metabolites and reactions, respectively. In particular, the red, 
blue and yellow nodes represent the sources, the sinks and the 
intermediate metabolites of the network, respectively. Red 
edges and blue edges represent the strictly mass-irreversible 
reactions, and violet red edges and navy blue edges represent 
mass-irreversible reactions. In particular, reactions in red or 
violet red proceed in the same direction as the original def- 
inition of the network [? ]; reactions in blue and navy blue 
proceed in a reversed direction. An infeasible reaction cycle 
is shown at the right top corner of Fig. 4. For simplicity, the 
currency of the biochemical networks, such as ATP, ADP, Pi, 
NADH, NAD, H+, H 2 and CO2, have not been shown. 



[? ] , but all other reactions are assumed to be reversible 
in default. In this network model, three transport fluxes 
are constrained to operate in the direction of net up- 
take of substrates (glucose, oleate and oxygen) and 28 
transport fluxes (including citrate, acetate and 26 addi- 
tional output fluxes) are constrained to operate in the 
direction of net production of reactants. A total of 33 
irreversibility constraints are specified and the directions 
of the remaining 217 internal reactions are computed. 



Given finite uptakes of glucose, oleate and oxygen, 
mass-balance analysis identifies 18 sMI internal fluxes 
(see red and blue fluxes in Fig. 4) , which represent a set 
of reactions essential for this network. It also shows that 
only 14 products (of 32 desired metabolites) could be pos- 
sibly produced from the sources (i.e. glucose, oleate, and 
O2), with the aids of three reversible transport fluxes 
of CO2, H + , and H2O. In addition, 70 mass-infeasible 
fluxes, 73 mass-irreversible fluxes, and 74 mass-reversible 
fluxes are found. Furthermore, the results of intensively 
sampled T-feasible networks (559,518,687 samples) fur- 
ther indicate that among 74 mass-reversible fluxes, three 
fluxes (CITRtm, ORNt3m 7 ORNUm) are T-infeasible, 
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FIG. 5: Demonstration of reductions in the number of re- 
action cycles with extra irreversibility constraints (EI). It 
represents the results of 10 5 independent simulations for gen- 
erating T-feasible solution spaces. Bars represent the stan- 
dard deviation of the number of cycles that are eliminated at 
each particular assignment. 

because these three reactions 

CITRtm : citrL[m] ^ citrL[c], 
ORNtim : h[c] + orn[m] ^ h[m] + orn[c], 
ORNtAm : citrL[c] + h[c] + orn[m] ^ citrL[m] + h[m] 
+orn[c] (3) 

develop a reaction cycle as CITRtm + ORNHm — 
RNtim = 0. It is also found that h[c] is the only 
metabolite that connects this reaction cycle to the re- 
maining part of the network, which reminds of the reac- 
tion cycle J2 — ► J3 — > J4 — > J2 in the example network. 
Thus neither direction is feasible for all these three reac- 
tions. 

The solution space defined by MIUTI needs to be fur- 
ther refined since 1099 cycles are still embedded in this 
network. However, by specifying only 14 reversible fluxes 
(out of 71 T-reversible fluxes) be irreversible, it is possi- 
ble to eliminate all those cycles. A final T-feasible net- 
work for mammalian cardiomyocyte has been illustrated 
in Fig. 4, which consists of 112 infeasible reactions, 101 
irreversible reactions and 44 reversible reactions. How- 
ever, such a linearly constrained subspace is not unique, 
as indicated in the example network. To further explore 
the nonlinear T-fcasible space, a large number of sub- 
spaces (see Fig. 5) have been generated by assigning 
directions to those T-reversible fluxes. Importantly, each 
irreversibility assignment is generated in such a way that 
the maximum number of reaction cycles can be elimi- 
nated. Therefore, the generated set of linear constraints 
defines a feasible solution space which requires a minimal 
number of irreversibility constraints, and thus represents 
the largest linear subspace in Fig. 2. 

Surprisingly, the computed T-feasible irreversibility 
constraints show significant difference when compared to 
the originally assigned directions, by which the space de- 



TABLE I: Comparison of a T-feasible network with high flexi- 
bility to the original network as well as its mass-balanced and 
thermodynamically balanced networks. 





PI 


PI+FBA 


PI+FBA+EBA 


EBA" 


Infeasible fluxes 





149 


149 


112 


Positive 6 fluxes 


156 


60 


61 


42 


Negative 6 fluxes 


3 


21 


21 


59 


Reversible fluxes 


98 


27 


26 


44 



"A T-fcasible network corresponds to the network shown in Fig. 
4. 

6 The positive reaction direction is defined based on the lcft-to- 
right reaction direction in the original model definitions [? ] . 



fined is not thermodynamically feasible (7 reaction cy- 
cles are included). To make it physically feasible, three 
mass-reversible fluxes (CITRtm, ORNtSm, ORNtAm) 
have to be eliminated, and at least one additional di- 
rection assignment, either specifying AKGMALtm (pre- 
viously assigned as reversible) be positive or specify- 
ing ASPTAm be negative, has to be added into the 
original definition of the network. Although these ad- 
ditional constraints do not invoke significant impacts 
on the overall results from [? ], they do change the 
flux distribution with the given boundary fluxes. While 
the revised network represents a special case of the T- 
feasible solution space, it is significantly different from 
the network which shows the maximum level of flexibil- 
ity (i.e. containing the maximum number of reversible 
reactions). Table 1 shows the difference between an ex- 
ample of the largest T-feasible subspaces (shown in Fig. 
4), and the original definition of the network (PI) and its 
mass-balanced network (PI+FBA) and thermodynami- 
cally balanced network (PI+FBA+EBA). The significant 
difference between the ab initio prediction and the pre- 
viously assigned directions, can be explained by the fact 
that mass-conservation and thermodynamic constraints 
are by far not the only constraints that a network must 
satisfy. Those non-identified may come from gene regula- 
tory constraints, topological constraints, environmental 
constraints [? ], and constraints that are accumulated 
from evolution and natural selection. By systematically 
studying the difference between the results from ab ini- 
tio predictions and experimental observations, we may 
find some interesting knowledge "gap", leading to more 
accurate predictions on the flux distributions. 

In summary, a biophysical constraint is applied to 
genome-scale biochemical networks for computing flux 
directions and determining a linearly constrained, mass- 
and thermodynamically-balanced solution space. Given 
a reaction network and a set of prior constraints on cer- 
tain flux directions (such as known transport flux direc- 
tions) the algorithm determines a minimal set of irre- 
versibility constraints (reactions directions) that is suffi- 
cient to guarantee T-feasibility. The developed methods 
could be used as a tool to systematically design a bio- 
chemical network, greatly facilitate the further flux bal- 
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ance analysis and energy balance analysis, improve their studied. Numerically, a linear solution space will make 
predictions of the flux distribution, and help capturing flux optimization more reliable and efficient, 
the real physically meaningful behavior of the network 



